function [results] = doppler_effect_analysis(varargin)
% 多普勒效应分析
% 功能：分析时变信道、时间相关性、多普勒功率谱、相干时间
% 输入参数：
%   - time_duration: 仿真时间 (秒), 默认 1
%   - sampling_rate: 采样频率 (Hz), 默认 1000
%   - speeds_kmh: 移动速度数组 (km/h), 默认 [10, 50, 100]
%   - carrier_freq: 载波频率 (Hz), 默认 2.4e9
%   - Nt: 发送天线数量, 默认 4
%   - Nr: 接收天线数量, 默认 4
% 输出：
%   - results: 结构体，包含多普勒效应分析结果
%
% 作者：周勇
% 日期：2024年

%% 参数解析
p = inputParser;
addParameter(p, 'time_duration', 1, @isnumeric);
addParameter(p, 'sampling_rate', 1000, @isnumeric);
addParameter(p, 'speeds_kmh', [10, 50, 100], @isnumeric);
addParameter(p, 'carrier_freq', 2.4e9, @isnumeric);
addParameter(p, 'Nt', 4, @isnumeric);
addParameter(p, 'Nr', 4, @isnumeric);
parse(p, varargin{:});

params = p.Results;

%% 添加路径
addpath('../Common');

%% 获取颜色定义
colors = color_definitions();

%% 基本参数计算
time_axis = 0:1/params.sampling_rate:params.time_duration;
speeds_ms = params.speeds_kmh * 1000 / 3600;
wavelength = 3e8 / params.carrier_freq;
doppler_freqs = speeds_ms / wavelength;

fprintf('=== 多普勒效应分析 ===\n');
fprintf('仿真时间: %.1f 秒\n', params.time_duration);
fprintf('采样频率: %d Hz\n', params.sampling_rate);
fprintf('移动速度: %s km/h\n', mat2str(params.speeds_kmh));
fprintf('载波频率: %.1f GHz\n', params.carrier_freq/1e9);

%% 时变信道分析
temporal_correlation = zeros(length(time_axis), length(params.speeds_kmh));
theoretical_corr = zeros(length(time_axis), length(params.speeds_kmh));

for speed_idx = 1:length(params.speeds_kmh)
    speed = speeds_ms(speed_idx);
    doppler = doppler_freqs(speed_idx);
    fprintf('\n分析速度 %d km/h (多普勒 %.1f Hz)...\n', ...
        params.speeds_kmh(speed_idx), doppler);
    
    % 生成时变信道
    H_time_varying = generate_time_varying_channel(params.Nt, params.Nr, time_axis, speed, params.carrier_freq);
    
    % 计算时间相关性
    H0 = H_time_varying(:, :, 1); % 初始信道
    for t = 1:length(time_axis)
        Ht = H_time_varying(:, :, t);
        % 计算归一化相关性
        correlation = abs(sum(sum(conj(H0) .* Ht))) / (norm(H0, 'fro') * norm(Ht, 'fro'));
        temporal_correlation(t, speed_idx) = correlation;
        
        % 理论相关性 (Jakes模型)
        theoretical_corr(t, speed_idx) = besselj(0, 2 * pi * doppler * time_axis(t));
    end
end

%% 多普勒功率谱计算
freq_axis = -200:1:200; % Hz
psd_matrix = zeros(length(freq_axis), length(params.speeds_kmh));

for speed_idx = 1:length(params.speeds_kmh)
    doppler = doppler_freqs(speed_idx);
    % Jakes谱
    psd = zeros(size(freq_axis));
    for f = 1:length(freq_axis)
        if abs(freq_axis(f)) <= doppler
            psd(f) = 1 ./ (pi * doppler * sqrt(1 - (freq_axis(f)/doppler)^2));
        end
    end
    psd_matrix(:, speed_idx) = psd;
end

%% 相干时间计算
coherence_times = zeros(1, length(params.speeds_kmh));
for speed_idx = 1:length(params.speeds_kmh)
    % 定义相关性下降到0.5的时间为相干时间
    coherence_idx = find(temporal_correlation(:, speed_idx) < 0.5, 1);
    if ~isempty(coherence_idx)
        coherence_times(speed_idx) = time_axis(coherence_idx);
    else
        coherence_times(speed_idx) = NaN;
    end
end

%% 可视化结果
figure('Name', '多普勒效应分析', 'Position', [250, 250, 1200, 800]);

% 时间相关性
subplot(2,2,1);
for speed_idx = 1:length(params.speeds_kmh)
    plot(time_axis, temporal_correlation(:, speed_idx), ...
         ['-', colors(speed_idx)], 'LineWidth', 2);
    hold on;
end
grid on;
xlabel('时间 (秒)');
ylabel('时间相关性');
title('信道时间相关性');
legend(arrayfun(@(x) sprintf('%d km/h', x), params.speeds_kmh, 'UniformOutput', false));

% 理论相关性 (Jakes模型)
subplot(2,2,2);
for speed_idx = 1:length(params.speeds_kmh)
    plot(time_axis, theoretical_corr(:, speed_idx), ...
         ['--', colors(speed_idx)], 'LineWidth', 2);
    hold on;
end
grid on;
xlabel('时间 (秒)');
ylabel('理论相关性');
title('Jakes模型理论相关性');
legend(arrayfun(@(x) sprintf('%d km/h', x), params.speeds_kmh, 'UniformOutput', false));

% 多普勒功率谱
subplot(2,2,3);
for speed_idx = 1:length(params.speeds_kmh)
    plot(freq_axis, psd_matrix(:, speed_idx), ...
         ['-', colors(speed_idx)], 'LineWidth', 2);
    hold on;
end
grid on;
xlabel('频率 (Hz)');
ylabel('功率谱密度');
title('多普勒功率谱');
legend(arrayfun(@(x) sprintf('%d km/h', x), params.speeds_kmh, 'UniformOutput', false));

% 相干时间分析
subplot(2,2,4);
bar(coherence_times);
grid on;
title('相干时间 vs 移动速度');
xlabel('速度配置');
ylabel('相干时间 (秒)');
set(gca, 'XTickLabel', arrayfun(@(x) sprintf('%d km/h', x), params.speeds_kmh, 'UniformOutput', false));

%% 结果打包
results = struct();
results.temporal_correlation = temporal_correlation;
results.theoretical_corr = theoretical_corr;
results.psd_matrix = psd_matrix;
results.coherence_times = coherence_times;
results.time_axis = time_axis;
results.freq_axis = freq_axis;
results.params = params;

fprintf('\n多普勒效应分析完成！\n');

end

%% 辅助函数
function H_time = generate_time_varying_channel(Nt, Nr, time_axis, speed, carrier_freq)
% 生成时变MIMO信道
wavelength = 3e8 / carrier_freq;
max_doppler = speed / wavelength;

num_time_samples = length(time_axis);
H_time = zeros(Nr, Nt, num_time_samples);

% 初始信道
H_time(:, :, 1) = sqrt(0.5) * (randn(Nr, Nt) + 1i * randn(Nr, Nt));

% 时变信道 (基于Jakes模型)
for t = 2:num_time_samples
    time_diff = time_axis(t) - time_axis(t-1);
    
    % Jakes相关性
    correlation = besselj(0, 2 * pi * max_doppler * time_diff);
    
    % 生成相关的新信道
    innovation = sqrt(0.5) * (randn(Nr, Nt) + 1i * randn(Nr, Nt));
    H_time(:, :, t) = correlation * H_time(:, :, t-1) + sqrt(1 - correlation^2) * innovation;
end
end